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Abstract 

2^ \ The paper combines theoretical and appUed ideas which have been previously 

considered separately into a single set of evolution equations for Numeri- 
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Vh ' cal Relativity. New numerical ingredients are presented which avoid gauge 

pathologies and allow one to perform robust 3D calculations. The potential 
of the resulting numerical code is demonstrated by using the Schwarzschild 
f^ , black hole as a test-bed. Its evolution can be followed up to times greater 

^vq ! than one hundred black hole masses. 
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I. INTRODUCTION 

In spite of the valuable work made by many Numerical Relativity groups all around the 
world, three-dimensional (3D) gravitational systems are still challenging (see Ref |jl| for a 
recent review of the field) . This is not only because the present generation of computers 
does not allow yet to perform high resolution 3D calculations. This is also, or at least we 
believe it is, because there are new numerical problems associated with 3D gravitational 
calculations which were not recognized in the previous work with spherically symmetric 
(ID) or even axially symmetric (2D) systems. To identify these new problems and to deal 
with them requires new developments, both from the theoretical and from the applied point 
of view. 

The theoretical developments, as presented in Section II, are made in the framework of 
the well known 3+1 formalism |0-^. We propose first to change the standard set of evolution 
equations (obtained from the space components of the four-dimensional Ricci tensor) to the 
one obtained from the space components of the Einstein tensor (which we called 'Einstein 
system' |^ for obvious reasons). 

We will now incorporate to this second-order system the idea of transforming the mo- 
mentum constraint into an evolution equation for a set of three supplementary quantities, 
which we used in former works with first-order evolution systems [§|,0]. Although the origi- 
nal idea was devised to obtain hyperbolic evolution systems, variants of the same approach 
have been considered recently p-pHj] as a tool to improve the results in specific numerical 
simulations. 

These modifications do change completely the causal structure of the original system, as 
discussed in Refs [§,0. The ones presented in Section III, on the contrary, are minor changes 
from the theoretical point of view but they seem to be crucial from the computational one. 



They address the problem of the so called 'gauge pathologies' [|Tl|,|T^], which consists in the 
explosive behaviour of the quantities associated with the gauge modes. This problem seldom 
happened in ID calculations (at least for the standard settings) but, in our experience, it 
turns out to be the rule in 3D calculations. 

A possible remedy to this problem could be trying to fix the gauge-related degrees of 



freedom by using the well known maximal slicing [|T3|. This requires, however, to solve a 3D 
elliptic equation at every time step, which implies a huge computational cost, specially in 
high resolution calculations. But this is not the only alternative available: the gauge-related 
degrees of freedom can be isolated from the rest of the system by using the well known 
conformal decomposition of the space metric [0,|14| as in the recent work of Shibata and 
Nakamura |Q , which use the harmonic slicing ||r5| , p!6|l , which is much less expensive from the 
computational point of view. 

Both alternatives, however, are not sufficient by themselves to allow for long term 3D 
calculations in strong field situations, like black hole space-times. Even if one tries to combine 
both approaches, the errors in gauge modes show a tendency to grow without bound which 
should be tackled using improved tools \T^. We have devised numerical counter-terms which 



have shown to be very effective in controlling these errors without modifying the analytical 
system of equations. As they are not restricted to maximal slicing, one is then allowed to 
use other slicings with a great increase in the code speed. The use of these counter-terms 
opens then the door to performing robust 3D numerical calculations at a level similar to 



that of previous 2D or even ID results. 

We present in Section IV our results for the single black hole case as a test-bed for 
our finite difference code. We obtain long term evolution simulations, up to more than 
one hundred times the black hole mass. A convergence analysis shows that the collapse 
speed one obtains (the rate at which the numerical grid falls into the hole) depends on the 
grid resolution, but it converges to the well established results for spherically symmetric 
(ID) black holes. The same is true for the preservation of the black hole mass outside the 
horizon. A finer analysis, however, shows that the outer boundary conditions do not preserve 
the constraints and this affects the consistency of the long term results. In spite of that, the 
code runs forever (although we obviously stopped it when all the numerical grid fell into the 
hole). 

II. SPACE PLUS TIME DECOMPOSITION 



We will use the well known 3+1 decomposition IQ-Q of the space-time metric, namely 

ds"^ = -a^ dt^ + 7,j {dx' + (3' dt) {dx^ + (3^ dt) , (1) 

where the lapse function a and the shift /3* are related to the choice of the space-time 
coordinates and 7ij is the induced metric on every time slice. The extrinsic curvature Kij 
(second fundamental form) of the slices is given by 

{dt-Cp)-ii, = -2aKi,, (2) 

where C stands for the three-dimensional Lie derivative operator. 

Four of the ten Einstein field equations do not contain time derivatives when expressed 
in terms of the previously defined quantities. These four constraint equations can be easily 
identified: the "energy", or Hamiltonian constraint 

2^2 G°° = ^'^R - tr{K^) + {tr Kf , (3) 

where ^^^R\s the trace of the three-dimensional Ricci tensor, and the "momentum constraint" 

aG\ = VkK\-d,{trK), (4) 

where index contractions and covariant derivatives are with respect to the induced metric 

lir 

The time evolution of Kij is given by a system obtained from the remaining Einstein 

field equations. For instance, the space components of the four-dimensional Einstein tensor 

can be written 

{dt - Cfs) Kij = -Wiaj + 

a f^Rij - 2Kfj + trK Kij - dj] 

-a/A -fij f^R - tr{K^) + {tr Kf - 2 tr G] , (5) 

where '•^^-Rjj stands for the three-dimensional Ricci tensor. The matter terms in the perfect 
fluid case can be computed from 



Gij = Svr [(/i + p)uiUj +pjij] , (6) 

where fi is the total energy density of the fluid, p is the pressure, and m, is its fluid 3- velocity. 
We call this set of equations (0,|) the Einstein evolution system. Notice that it differs 
from the standard 3+1 evolution system (obtained from the space components of the four- 
dimensional Ricci tensor) by a term containing the Hamiltonian constraint. This ensures 
that the matter contribution (|^) vanishes in the Newtonian limit. Even in the vacuum 
case, the Einstein evolution system has proved to be superior when following the long term 
evolution for a spherically symmetric (ID) black hole 0. 



A. The momentum constraint as an evolution equation 

Let us consider now the quantities Vi, which can be obtained from the metric first 
derivatives as follows 

V, = l/2Y'{d^irs-dri^s). (7) 

One can alternatively compute the time evolution of Vi by using the momentum constraint 
(I) as a dynamical equation, that is ^■. 

dtV, - (3'dkV - (0^(3^) Vr + 

a^ G\ + ar K\ -UitrK- 

a [2 K[ Vr + K"-' r„,, - K[ dr{ln^]. (8) 

We have chosen this second approach, which amounts to consider the following set of basic 
dynamical quantities 

{7,,, Ki,, V}, (9) 

so that the condition (^ can be now considered just as an algebraic constraint which will 
hold if and only if the momentum constraint is satisfied. 

The three-dimensional Ricci tensor can be written in terms of these quantities as 

(3)i?,,. = -1/2 7'-^ dl 7i, + 5,5 In^ - d, V, - d, V + 
rf^. (2 Vk - dk In^) + 7'' 7^' m lri){di 7,,) 

— Tjfcr Tjls] , (10) 

so that its trace can be easily expressed as follows 

(3)/? = -2 dk V^ + Tkrs r^^^^ - 7''^ {dr In^) (9, In^) . (11) 



B. Gauge conditions 

So far we have considered a completely general coordinate system. We will keep our 
freedom to fix the shift components /5* as purely kinematical quantities. The choice of the 
lapse function a, on the contrary, will be closely related with the dynamics of our system. 

A popular choice is to impose the maximal slicing condition |T3| 



trK = 0, (12) 

which amounts to compute a by solving the elliptic equation obtained by taking the trace 
of d^) or a combination of this trace with the Hamiltonian constraint @. This choice has 
proven to be quite effective in some cases, although it is not free from long term instabilities. 
Also, elliptic equations are very expensive in terms of computational resources and getting 
high resolution long term calculations on three-dimensional grids seems to be beyond the 
capabilities of present day computers. 

A different strategy is to link the evolution of the lapse a to that of the volume element 
^/7 by means of the following equation 

{dt — Cfs) Ina = —a f tr K, (13) 

where / is a given non-negative function of a. This opens the way to many different choices 
of the time slicing. The geodesic slicing is included as a subcase with / = 0. The / = 1 case 
corresponds to the "harmonic slicing" |T5|JT6| (the resulting time coordinate is harmonic). 
Another interesting case is obtained when / = 1/a; it mimics maximal slicing near a singu- 



larity, when the lapse collapses to zero (a very similar case was considered in |T^|T^. The 
related choices / = n/a {n = 2, 4) mimic both the singularity- avoidance and the large- 
gradients-avoidance properties of the maximal slicing for a fraction of the computational 
cost. 

So far we have obtained a closed dynamical system (once the shift /5* is given) which 
can be considered as the second-order equivalent of the first-order hyperbolic system given 
in Ref 0,11. The causal structure of this second-order system will be discussed elsewhere 
p^ . All we can say by now is that the direct equivalence with a first-order system which 
is known to be hyperbolic seems a good starting point for numerical relativity applications. 
But this feeling should be confirmed by experience, as we will see in what follows. 

III. A ROBUST EVOLUTION SYSTEM 
A. Keeping apart trace and trace-free components 



Both from the theoretical and from the practical point of view PJT^, it is clear that 



tr K is directly related to the gauge degrees of freedom, which should be treated very 
carefully to avoid spurious 'gauge pathologies' or riddles that do appear too frequently 



in 3D calculations [|I2[. This suggests the convenience of evolving separately the trace and 
trace- free components of Kij. 

We will follow here the notation of Shibata and Nakamura f^ to write down the following 
conformal decomposition of the space metric 



4$ ~ 



(14) 



with $ chosen so that the determinant 7 of the conformal metric is equal to unity, namely 

V7 = e^* . (15) 

This means that we can split out the evolution equation (^ into the following pieces pHtJ : 

{dt - Cp) %j = -2a Aij , (16) 



{dt~Cf3)^ = -1/6 a K 



where we have noted for short 

K = trK. 
We can also split out the evolution equations (^, H) in the same way: 

{dt - Cf,)Aij = e-^^i-Via, + 2 {a^ $, + aj ^i) 

+ l/3r'(V.as-4«,$,)7ii 

+ a l^^^Ri, - 1/3 {fr ^^^^R) 7^,- - dj ] } 

+ a{K Aij-2f''AirAjs) 



(17) 



(18) 
(19) 



(20) 



{dt - Cp)K = e-4*{-(7-V.a.) - 2 (^^a, $,) 
+ a [1/4: fr^^^R+ 1/2 frG] } 
+ a{l/2K^ + ?,/Ar'f'AkrAi, 



Us 



(21) 



dtV, - p'^dkV, - {d,Pn Vr + 1/2 (7'^7.s - <5f S:) dlP' 
= «' G° + AirY'{as -2aVs + 2a $,) 
-a7"^7''^fc;f,,,-2/3ira, . 



(22) 



The resulting system (0, |T^, ^, ^, ^) allows to compute the time evolution of the 
enlarged set of dynamical quantities 



{-fij , $ , Aij , K , Vi} 



(23) 



which has the same basic structure of the Einstein evolution system (|, |^, ||) because the 
two extra quantities we have introduced are redundant, so that the constraints 



detj = l, trA = Q 



(24) 



are trivially propagated in time. 



B. Numerical counter-terms 



In numerical applications the algebraic constraints ( P^ can be easily monitored as error 
indicators. Their behaviour is quite unstable: tr A grows without bound (the actual rate 
depends on every case) and this makes det 7 to go far away of its exact value of 1. This leads 
to an obvious inconsistency which causes the numerical 3D code to crash after a relatively 
short time. 

To fight against these errors, we have modified the right-hand-sides of the evolution 
equations (|16| , ^0]) by adding the following counter-terms: 

(dt - Cp) 7,, = ... - 1/3 {det 7 - l)/r 7,,- , (25) 

{dt - Cp) \^= ... - 1/3 {tr A) It % , (26) 

where r is a 'time constant'. 

It is easy to see that the counter-terms (|25| , |26| ) by themselves would cause an exponential 
decrease of the errors in tr A and det 7. In practice, we have set the time constant r to the 
same value of the numerical time step dt in order to keep these errors within an small order 
of magnitude. 

One should notice at this point that these counter-terms do vanish in the continuum 
limit (at least for second-order-accurate numerical algorithms), so that they do not affect 
at all the structure of the analytical evolution system. Harmless from the theoretical point 
of view, they provide a valuable tool that allows us to perform robust 3D calculations, as 



we will see below. See Refs. |l21| , |22[| for the use of numerical control terms in elliptic gauge 
conditions. 



IV. APPLICATION TO THE BLACK HOLE CASE 

Let us start with the usual time symmetric initial data {Kij = 0) for a conformally flat 
metric (space isotropic coordinates), corresponding to a Schwarzschild black hole. Let us 
now replace the interior region by a constant density incoherent matter content, keeping the 
exterior vacuum region of the original hole |^. This 'stuffed' black hole approach allows 
one to start from singularity-free initial data and avoids then the usual procedure of excising 
the inner region from the computational grid, which would introduce an internal boundary 
at the black hole horizon. 

We have performed our calculations in the Cactus code []I| environment. In order to 
save computational resources, we have taken advantage of the equatorial symmetry of the 
problem to evolve just an octant of the space slices. This requires to introduce boundary 
conditions at the symmetry planes. The outer boundary has been treated as in Ref ^. In 
Fig. m we consider the effect of this boundary condition over the collapse of the lapse. The 
results seem to be quite independent on the position of the outer boundary. 



FIGURES 




FIG. 1. The evolution of a is compared for two positions of the external boundary (r = lOM 
and r = 20M) for a Schwarzschild black hole. The resolution is set at dx = 0.3M. The lapse 
is plotted at t = 21M, 37M and 47M. In the latter the black hole has gone almost out of the 
smallest grid. We can see that the agreement is quite good. 

As coordinate conditions, we use normal coordinates (/5* = 0) and an algebraic slicing 
of the form (|T^ with the choice / = 2/a which leads to smoother profiles of the evolved 
quantities. To start a dynamical simulation we can take as usual a constant initial lapse. The 
time evolution is then obtained by using the 'icn' second-order-accurate algorithm which has 
been incorporated to the Cactus version currently used by the Potsdam group. The actual 
plight for the reliability of this collapse turns out to be resolution, as we can see in Fig. |^. 
The 3D results converge to the well established ID solution and we see that a resolution 
finer than O.IM is needed if one wants to get really close to the physical solution. 




FIG. 2. The convergence of the lapse towards ID results is shown, a is plotted against r/M 
at t = 30M for a Schwarzschild black hole. The solid line corresponds to the ID result with 
high resolution. The remaining lines, converging to that solution, correspond to 3D results with 
dx = 0.3M, O.lSAf and 0.075M. 



As far as our test-bed problem has spherical symmetry, we can make use of the Bondi 
mass function to monitor the convergence of the numerical solution towards the analytical 
one. In Fig. ^ we plot the time evolution of the numerically obtained mass at the apparent 
horizon. The results show the expected convergence to the constant line corresponding to 
the analytical value. This confirms our previous conclusions. 




FIG. 3. The Bondi mass at the apparent horizon is plotted against time for a Schwarzschild 
black hole. Results with three resolutions (dx = 0.3M, 0.15M and 0.075M) are plotted, which 
show convergence to the analytical value itiah/M = 1. 

Finally, we show in Fig. | the performance of the code in long runs. We can follow the 
black hole evolution for unlimited time. We just stop the calculation when the numerical 
grid is almost completely inside the hole. 




FIG. 4. The collapse of the lapse for a Schwarzschild black hole is shown, a is plotted against 
the radius for increasing times (t = 15M, 30M, 45M, 60M, 75M, 90M, 105M, 120M). The whole 
grid goes eventually inside the hole and the calculation is stopped (although it never crashes). At 
late times, however, the numerical solution departs from the physical one due to boundary effects. 



V. CONCLUSIONS 

The equations we propose are interesting both from the theoretical and for the numerical 
point of view. The decomposition of Kij into trace and trace-free parts allows a separate 
monitoring of the numerical errors which leads to the following quite surprising conclusion: 
the numerical trace of the trace-free part becomes easily unstable in strong field 3D calcula- 
tions. Notice that this decomposition is not made in most Numerical Relativity formalisms, 
so that the errors arising from the trace-free part are directly mixed up with the dynamics 
of the trace part. This might explain the arising of the so called 'gauge pathologies' |1TT|JT^ 
even in situations with apparently harmless initial data. 

The counter-terms we have introduced allow one to deal with this problem and this 
is enough to perform robust 3D Numerical Relativity calculations, even in the black hole 
case. We have demonstrated this by actually running the code up to times greater that one 
hundred black hole masses without any sign of instabihty: we just stopped when nearly all 
the numerical grid was inside the hole. 
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